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Abstract 

We introduce a model for ionic electrodiffusion and osmotic water flow through 
cells and tissues. The model consists of a system of partial differential equa- 
tions for ionic concentration and fluid flow with interface conditions at deforming 
membrane boundaries. The model satisfies a natural energy equality, in which 
the sum of the entropic, elastic and electrostatic free energies are dissipated 
through viscous, electrodiffusive and osmotic flows. We discuss limiting models 
when certain dimensionless parameters are small. Finally, we develop a numer- 
ical scheme for the one-dimensional case and present some simple applications 
of our model to cell volume control. 



1. Introduction 



Systems in with important electrodiffusion and osmotic water flow are found 
throughout life [H H, H| . Such systems include brain ionic homeostasis 0, 5 1 , 
fluid secretion by epithelial systems [6| , electrolyte regulation in the kidney [7|, \8 1, 



ui. 



fluid circulation in ocular systems [9|, [f0[ , and water uptake by plants 

Mathematical models of electrodiffusion and/or osmosis have been proposed 
and used in many physiological contexts, and have formed a central topic in 
biology for a very long time [3, [lil [f3| . Some are simple models using ordinary 
differential equations while others are more detailed in that they include partial 
differential equations (PDEs) descri bing the spatial variation of the concentra- 
tion and flow fields [3 El, [3 El , 18, U HI- I n this paper, we propose a system 
of PDEs that describes ionic electrodiffusion and osmotic water flow at the cel- 
lular level. To the best of the authors' knowledge, this is the first model in which 
osmotic water flow and electrodiffusion have been treated within a unified frame- 
work including cells with deformablc and capacitance-carrying membranes. A 
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salient feature of our model is that it possesses a natural thermodynamic struc- 
ture; it satisfies a free energy equality. As such, the present work may be viewed 
as a generalization of the classical treatment of osmosis and electrodiffusion in 
irreversible thermodynamics to spatially extended systems [21, 22, 23 1. 

To introduce our approach, wc first focus attention on uncharged systems. 
In Section [3J we treat the case in which the diffusing chemical species carry no 
electric charge. We write down equations that are satisfied by the water velocity 
field u, the chemical concentrations Cfe, k = 1, • • • , N and the membrane position 
X. The model is shown to satisfy a free energy equality in which the sum of the 
entropic free energy and the elastic energy of the membrane is dissipated through 
viscous water flow, bulk diffusion, transmembrane chemical fluxes and osmotic 
water flow. One interesting consequence of this analysis is that the classical 
van t'Hoff law of osmotic pressure arises naturally from the requirement that 
osmotic water flow be dissipative. We note that models with the similar purpose 
of describing diffusing non-electrolytes and their interaction with osmotic water 
flow across moving membranes, have been proposed in the literature 1||, 24, 25| 
(in the App endix |Appendix A. 2 wc discuss the relationship of our model with 
that of [3 [Hi). 

In Section [3J we extend the model of Section [2] to treat the case of ionic 
electrodiffusion. We introduce the electrostatic potential 4> which satisfies the 
Poisson equation. The membrane now carries capacitance, which can result in 
a jump in the electrostatic potential across the membrane. We shall see that 
this model also satisfies a free energy equality. The free energy now includes 
an electrostatic contribution. The verification of the free energy equality in 
this case is not as straightforward as in the non-electrolyte case, and requires a 
careful examination of surface terms. 

In Section 01 wc discuss simplifications of our model. We make the system 
dimcnsionless and assess the relative magnitudes of the terms in the equations. 
An important simplification is obtained when we take the clcctroncutral limit. 
In this case, the electrostatic potential becomes a Lagrange multiplier that helps 
to enforce the electroneutrality condition. 

In Section O we develop a computational scheme to simulate the limiting 
system obtained in the electroneutral limit, when the geometry of the cell is 
assumed spherical. As an application, we treat animal cell volume control. 



2. Diffusion of Non-electrolytes and Osmotic Water Flow 

2.1. Model Formulation 

Consider a bounded domain C K 3 and a smooth closed surface T C 0. 
This closed surface divides fi into two domains. Let Vti C £1 be the region 
bounded by L, and let Vt e = f2\(f2jLir). In the context of cell biology, f2, ; may be 
identified with the intracellular space and Q e the extracellular space. Although 
cell physiological systems of biological cells serve as our primary motivation for 
formulating the models of this paper, this identification is not necessary. 
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In this section wc formulate a system of PDEs that governs the diffusion of 
non-electrolytes and osmotic flow of water in the presence of membranes. In 
Section [31 we shall build upon this model to treat the electrolyte case. 

We consider N non-electrolyte chemical species whose concentrations we call 
Cfc, k = 1, ■ • • , N. Let uj be the entropic part of the free energy per unit volume 
of this solution. The following expression for u) is the most standard choice: 

N 

W = 2_j k B Tc k hi Cfc. (2.1) 

fc=l 

This expression is valid when the ionic solution is sufficiently dilute and lead 
to linear diffusion of solute. Our calculations, however, do not depend on this 
choice of w. If the solution in question deviates significantly from ideality, other 
expressions for u> may be used in place of ujq. 

Given u, the chemical potential fj,h of the k— th chemical species is given as: 

du . . 

fi k = o-fc, cr fc = — - . (2.2) 
oc k 

We have introduced two symbols fj,k and crfc in anticipation of the discussion of 
the electrolyte case, where and crfc are different. For water, it is convenient 
to consider the water potential ip w , the free energy per unit volume, rather than 
the fi w , the chemical potential (free energy per molecule). The water potential 
and water chemical potential are thus related by the relation v w ip w = fj, w where 
v w is the volume of water per molecule. We define: 



ip w = <k w + p, ir w = \uj - y CfctJfc = [uj ~y~] Cfc-— (2.3) 




where p is the pressure. As we shall see, p will be determined in our model via 
the equations of fluid flow (Eq. (|2.8|0 . The entropic part of rp w (or the osmotic 
pressure), tt w , is given as the negative multiple of the (scmi)Lcgcndrc transform 
of uj with respect to all the ionic concentrations Cfc . This expression for osmotic 
pressure can be found, for example, in [26j |. As we shall see in the proof of 
Theorem 12.11 the above definition of tt w is forced upon us if we insist that our 
model satisfy a free energy dissipation principle. 

We begin by writing down the equations of ionic concentration dynamics. 
At any point in Qi or f2 e 

^ + V.(uc k ) = V.(c^V«) (2.4) 

where is the diffusion coefficient and u is the fluid velocity field. Ions thus 
diffuse down the chemical potential gradient and are advected with the local fluid 
velocity. We have assumed here that cross-diffusion (concentration gradient of 
one species driving the diffusion of another species) is negligible. 
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We must supplement these equations with boundary conditions. Most for- 
mulations of non-equilibrium thermodynamic processes seem to be confined ei- 



ther to the bulk or to the interface between two bulk phases [2JJ, |27J, [23[ . Here 
we must couple the equations in the bulk and with boundary conditions at the 
interface, which as a whole give us a consistent thermodynamic treatment of 
diffusion and osmosis. 

On the outer boundary r out = dfl, for simplicity, we impose no-flux bound- 
ary conditions. Let us now consider the interfacial boundary conditions on the 
membrane T. Since we want to account for osmotic water flow, the membrane T 
will deform in time. Sometimes, we shall use the notation T t to make this time 
dependence explicit. Let F re f be the resting or reference configuration of T. The 
membrane will then be a smooth deformation of this reference surface. We may 
take some (local) coordinate system 8 on r re f, which would serve as a material 
coordinate for T t ■ The trajectory of a point that corresponds to 9 = 9q is given 
by X(0 o ,t) 6 R 3 . For fixed t, X(-,i) gives us the shape of the membrane T t . 

Consider a point x = X(0, t) on the membrane. Let n be the outward unit 
normal on F at this point. The boundary conditions satisfied on the intracellular 
and extracellular faces of the membrane are given by: 

( D k \ <9X , , 

Ck ( u - JT^f^Vk I • n = c k — ■ n + j k + a k on or T e . (2.5) 

The expression "on r i e " indicates that the quantities are to be evaluated on 
the intracellular and extracellular faces of T respectively. The term j k is the 
passive chemical flux that passes through the membrane and is the active 
flux. Fluxes going from Q, to fi e is taken to be positive. Equation (|2.5[) is just 
a statement of conservation of ions at the moving membrane. It is easy to check 
that (|2.4p together with (|2.5[) implies conservation of each species. 

The flux jk is in general a function of concentrations of all chemical species 
on both sides of the membrane. The chemicals are usually carried by channels 
and transporters, and the functional form of jk describe the kinetic features of 
these carriers. As we shall see in Section POl jk can also be a function of the 
difference in water chemical potential across the membrane. 

The passive nature of the jk is expressed by the following inequality: 

\Hk\jk > 0, \p,k] = A t fe| r! - Mfc| re (2.6) 

where -| r . expresses evaluation of quantities at the intracellular and extra- 
cellular faces of the membrane T respectively. We shall see that this condi- 
tion is consistent with the free energy identity (|2.20[) . For any quantity w, 
[to] = tw| r . — u>| r will henceforth always denote the difference between w eval- 
uated on the intracellular and extracellular faces of T. A simple example of jk 
occurs when jk is a function only of [fik] and satisfies the following conditions: 

Jk=3k(Wl), j*(0) = 0, >0. (2.7) 
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It is easily checked that jk in this case satisfies (|2.6|) . We shall see concrete 
examples in Section [3] Condition (|2.6p is somewhat restrictive in the sense that 
it is possible to have a "passive" flux that does not satisfy (|2.6p when multiple 
species flow and interact. We shall discuss this briefly in Section |2~51 Condition 
(|2.6I) needs to be relaxed to describe systems in which different chemical species 
flow through one channel or (passive) transporter. Those systems usually couple 
fluxes of different chemical species. They often couple (unidirectional) influx 
and efflux of the same species (symporters and antiporters) [H, [l2|, HH, 0, Q . 
The active flux ak is typically due to ionic pump currents often driven by ATP 

man. 

We now discuss force balance. We shall treat the cytosol as a viscous fluid 
and the cell membrane as an elastic surface. The cell membrane itself is just a 
lipid bilayer, and cannot support a large mechanical load. The cell membrane is 
often mechanically reinforced by an underlying actin cortex and overlying system 
of connective tissue, and in the case of plant cells, by an overlying cell wall. If wc 
view these structures as part of the membrane, our treatment of the membrane 
being elastic may be a useful simplification. Wc could employ a more complete 
model of cell mechanics incorporating, in particular, the mechanical properties 
of the cytoskeleton and extracellular lamina. However, our emphasis here is on 
demonstrating how osmosis can be seamlessly combined with mechanics, and we 
intentionally keep the mechanical model simple to clarify the underlying ideas. 

Consider the equations of fluid flow. The flow field u satisfies the Stokes 
equation at any point in fii or fi e : 

i/Au- Vp = 0, Vu = (2.8) 

where v is the viscosity of the electrolyte solution. Note that the above equations 
can also be written as follows: 

V-£ m (u,p) =0, V-u=0, S m (u,p) = j/(Vu+(Vu) T )-p/ (2.9) 

where I is the 3x3 identity matrix and (Vu) T is the transpose of Vu. Here, 
S m is the mechanical stress tensor. It is possible to carry out much of the 
calculations to follow even if we retain inertial terms and work with the Navier- 
Stokes equations or use other constitutive relations for the mechanical stress. 
In particular, such modifications will not destroy the free energy identity to be 
discussed below. We do note, however, that incompressibility is important for 
our computations. 

Wc now turn to boundary conditions. We let u — on the outer boundary 
Tout for simplicity. On the cell membrane T, we have the following conditions. 
Take a point x = X(0,t) on the boundary T, and let n be the unit outward 
normal on T at this point. First, by force balance, we have: 

[E m (u,p)] =F e i as . (2.10) 

Here, F c i as is the elastic force per unit area of membrane. 

We make some assumptions about the form of the elastic force. We assume 
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that the membrane is a hyperelastic material in the sense that the elastic force 
can be derived from an elastic energy functional -E c ias that is a function only of 
the configuration X: 

£ clas (X) = / E (X)dmr r8f (2.11) 
•/r ref 

where mr roI is the surface measure of r re f and £ is the elastic energy density 
measured with respect to this measure. It is possible that £ is a function of 
spatial derivatives of X. The elastic force F c i as (x) satisfies the relation: 



d_ 

ds 



f £(X(9) + sY(6))dmr tei = - [ F clas (x) • Y(X" 1 (x 
=o Jr TBt Jr 



))dm T (2.12) 



where Y is an arbitrary vector field defined on r ro f and mr is the natural 
measure on the surface T and is related to mr ref by dmr = Qdmr, a{ where Q 
is the Jacobian determinant relating Tt to the reference configuration r rc f. The 
expression X _1 (x) is the inverse of the map x = X(0). Thus, F e i as is given 
as the variational derivative of the elastic energy up to the Jacobian factor Q. 
Consequently, we have the following relation: 

-£elas(X) =~ J Folas ■ -Q^dm T . (2.13) 

In the above, ^ should be thought of as a function of x, i.e., ^ = ^r(X -1 (x)). 
We shall henceforth abuse notation and let ^ be a function of x or depending 
on the context of the expression. 

In addition to the force balance condition (|2.10[) . we need a continuity con- 
dition on the interface T. Since we are allowing for osmotic water flow, we have 
a slip between the movement of the membrane and the flow field. At a point 
x = X(0, t) on the boundary T we have: 

u-—=j w n (2.14) 

where j w is water flux through the membrane. We are thus assuming that water 
flow is always normal to the membrane and that there is no slip between the 
fluid and the membrane in the direction tangent to the membrane. Given that 
n is the outward normal, j w is positive when water is flowing out of the cell. 
Like jk in (|2.5[) . we let j w be a passive flux in the following sense. We let: 



jw = jw(bPw}), = k w \ t - ((£ m (u,p)) • n)| r . (2.15) 



where n w was the entropic contribution to the water potential defined in (|2.3|) . 
In general, j. w can be a function of other variables (see Section l2~3"|) . The function 
j w satisfies the condition, analogous to 



[ip w ]j w > 0. (2.16) 



G 



This condition is clearly satisfied if: 



jwibPw}), jw(o) = o, 



dj w 



> 0. 



(2.17) 



Jw — 



Water flow across the membrane is thus driven by the difference in the cntropic 
contribution to the chemical potential as well as the jump in the mechanical 
force across the cell membrane. When the flow field u is equal to 0, the above 
expression for tjj w reduces to: 



We may thus view ip w as a modification of i\) w to take into account dynamic 
flow effects. Let uj = luq given in (|2.1[) in the definition (|2.3|) of tt w . Then, we 
have, under zero flow conditions, 



We thus reproduce the standard statement that water flow across the membrane 
is driven by the difference in osmotic and mechanical pressure, where the osmotic 
pressure, tt w , is given by the van t'Hoff law. 

2.2. Free Energy Identity 

We now show that the system described above satisfies the following free 
energy identity. 

Theorem 2.1. Suppose Ck,u,p be smooth functions that satisfy (|2.4|) . (|2.8[) . 
and in VLi and VL e and satisfy boundary conditions (|2.5[) . (|2.10|) . (|2.14|) on the 

membrane V . Suppose further that Ck and satisfy no- flux boundary conditions 
and u = on the outer boundary T out . Then, Ck, u,p and (j) satisfy the following 
free energy identity. 

"77 (^S + E e i as ) = —I p — J p — J a 



a- 



TTu. | r . e + P\ r . e = tpw\ r 



(2.18) 




(2.19) 




(2.20) 
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Here, E e i as was given in (|2.11[) and |Vu| is the Frobenius norm of the 3x3 rate 
of deformation matrix Vu. 

// afe = 0, then the free energy is monotone decreasing. 

The free energy is given as a sum of the entropic contribution Gs and the 
membrane elasticity term £7 e las- This free energy is dissipated through bulk 
currents I p and membrane currents J p . Dissipation in the bulk comes from ionic 
clcctrodiffusion and viscous dissipation. Dissipation at the membrane comes 
from ionic channel currents and transmembrane water flow. If active membrane 
currents are present, they may contribute to an increase in the free energy 
through the term J a . It is sometimes useful to rewrite J p + I p as: 

Jp ~\~ Jp F w -\- F c , 

F w = / v\Vu\ 2 dyL+ j [tpwjjwdmr, 

Jn.un. Jr (2.21) 



f N D k 2 f N 

Fc.= y2 c k-r~ \VfJ-k\ 2 dx+ / yVfcbfcdmr, 



where F w and F c are the dissipations due to water flow and solute diffusion 
respectively. In the statement of the Theorem, it is important that (|2.6[) and 
(|2.16|) are used only to conclude that J p be positive. Identity (|2.20|) should be 
seen as giving us the definition of what a passive current should be. 

We now prove Theorem 13.11 An interesting point about the calculation 
to follow is how dissipation through transmembrane water flow comes from 
two different sources, equations for ionic concentration dynamics and the fluid 
equations. The former contributes the osmotic term tt w term, and the latter 
contributes the mechanical term p, which together add up to the water potential 

Proof of Theorem \2.1i First, multiply (|2.4[) with fik in (|3.1|) and integrate over 
fij and sum in k: 

Y,J n ^ (j£ + V- (uc k )j d x = J q MfeV • (cfc^V^j dx. (2.22) 
The summand in the right hand side becomes: 



MfcV ■ I Cfey^V^fc J dx= [ ( ^fcCfeT^V^fc ■ n I dm r 



n 4 \ k B T ) J Ti \ k B T 



/(^IV^)*. 



(2.23) 
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where n is the outward normal on T. Consider the left hand side of (|2.22[) . 



k=l 



( £ + V ■ W ) =E^-(t5T + u V «) = S+V-(-), (2.24) 



9w / 9cfc 



fe=i 



9cfe V at 



9a; 



dt 



where we used (|2.2[) and the incompressibility condition in (|2.8[) . Integrating 
the above over fli, we have: 



dt 



— + v • M Ux = / — dx 



9w 



at 



d 
dt 



wu • ndmr 

ax 
~a7 



£J U — 



(2.25) 



where we used the fact that u is divergence free in the first equality. The term 
involving ^ comes from the fact that the membrane T is moving in time. 
Performing similar calculations on f2 e , and adding this to the above, we find: 



— f u;dx+ I [uj]j w dm r 
dt Jniun,, Jr 



N 

E 

k=l 



MfcCfe— -V/i fc • n 
KrT 



N 

dmr — y I 

k=1 J(i, 



(2.26) 



where we used ()2.14j) . Using f|2.5|) and ()2.14|) . we may rewrite the second bound- 
ary integral as follows: 



L 



MfcC fe — — V/ifc ■ n 
k B T 



dm r = J ([jUfeCfc] j w - \jJ,k](jk + a k ))dm v . (2.27) 
We now turn to equation ()2.8|) . Multiply this by u and integrate over fi;: 
/ u • (uAu - Vp)dx = / (£ m (u,p)n) • udm r - / ^|Vu| 2 dx = (2.28) 

Jtti JTi JUi 

Performing a similar calculation on VL e and adding this to the above, we have: 



[(E m (u,p)n)] -udmr - / j/|Vu|-dx = (2.29) 

We may use (f27TO]) . (f27T3|) and ([2~T4| to find 

-^E clas (X)= f [{Z m (u,p)n)-n]j w dmr- I ^|Vu| 2 dx (2.30) 
at Jr JQiun s 



Combining ([2350 . (j2~2T)) and (f23Djl . we have: 



d 
dt 



/ wdx + £' e i as (X) 



TV 

E 



fe=i 



/ (c k ^-\Vfi k \ 2 )d^- f v\\7u\ 2 dx (2.31) 
- J \pk\(jk + a k )dm T — J [u- CfcfTfc - (E TO (u,p)n) • n] j w dm r 

Recalling the definition of ip w in (|2.15|) . we obtain the desired equality. In 
the absence of active currents a*, is decreasing given that jk and j w satisfy 
conditions (|2.6[) and f|2 . respectively. □ 

In the last line of the above proof, note that the expression for ip w arises 
naturally as a result of integrating by parts. In this sense, we may say that 
osmotic water flow arises as a natural consequence of requiring that the free 
energy be decreasing in time. 

2.3. Cross Coefficients and Solvent Drag 

As can be seen from (|2.20[) or (|3 . 1 5[) the only condition we need to impose 
for the free energy to decrease with time in the absence of active currents is the 
following: 

N 

[^»}ju, + J2^k]jk>0. (2.32) 

fc=i 

This condition is weaker than conditions (|2.6p and (|2.16[) being satisfied sep- 
arately by jk and j w . We now discuss an important case in which jk and j w 
may not individually satisfy (|2.6[) and (|2 . 16[) but (|2 .32[) is satisfied nevertheless. 
This arises whenever fluxes are coupled as is usually the case for fluxes through 



transporters or (single filing) channels [28|, [12|, [29|, |2j, [3 . We note that such 



31, 32, 33, 3 



cross-diffusion can be relevant even in bulk solution 

If \pk] and [ipw] remain small, the dissipation J p in (|2.20[) may be approxi- 
mated by a quadratic form in the jumps: 

J P = M • jdm r = / [fi] ■ (C [n])dm r , 

Jr Jr (2.33) 

V = (Ml, ■ • ■ ,VN,^w) T ,j = (jl, ■ ■ • ,jN,jw) T , 

where £ is a symmetric (N + 1) x (N + 1) matrix. Requiring that the free 
energy be decreasing implies that C must be positive definite. The maximum 
dissipation principle requires that j be given as variational derivatives of J p /2 
with respect to [/it]: 

j = L [fj] . (2.34) 
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Note that, without the maximum dissipation principle, (|2.33p only implies j = 
(C+C) [fi] where C is an arbitrary skew symmetric matrix. The symmetry of the 
coefficient matrix C relating [fi] and j is an instance of the Onsagcr reciprocity 
relation [H ED, HI | . 

A lipid bilayer membrane is impermeable to many solutes, but only approx- 
imately so. In this case, a water flux may induce a solute flux, and this may be 
expressed as Ck w ^ where Ckw is the (k, N+ 1) entry of the matrix C. This is 
known as solvent drag. Given the presence of such cross coefficients, (|2.6[) and 
(|2.16p are not necessary true, whereas condition (|2.32|) is true by construction. 



3. Electrodiffusion of Ions and Osmotic Water Flow 

3.1. Model Formulation 

Let us now consider the case in which the chemical species are electrically 
charged. As in the previous section, we let Ck, k = 1, • ■ ■ , N be the concentra- 
tions of the ionic species. Given ui, the entropic part of the free energy per unit 
volume, the chemical potential fj,k of the k— th species of ion is given as: 

fj>k = + Q z k4> = °~k + qzk4>- (3.1) 

The chemical potential is thus a sum of the entropic term o~k and the electrostatic 
term. In the electrostatic term, q is the elementary charge, is the valence of 
the A:-th species of ion, and <j) is the electrostatic potential. The definition of 
the water potential, ip w and ipw, remain the same. The ionic concentrations Ck 
satisfy (|2.4[) and (|2.5|) except that we now use (|3.1j) as our expression for the 
chemical potential. Ions are thus subject to drift by the electric field in addition 
to diffusion and advection by the local flow field. 

If the electrolyte solution is sufficiently dilute, the chemical potential /ife is 
given by (|3.1[) with ui equal to (|2.ip . However, deviations from ideality can be 



significant in electrolyte solutions, especially in higher concentrations [35l 1361 
[37L [3sL I39I ] . Cross-diffusion (or flux coupling) in the bulk can also be significant 
in electrolyte solutions [13, El, [H, [33 , 34 1 . These effects are clearly important in 



describing the molecular physiology of ion channel pores and enzyme active sites 
at which ionic concentrations can reach tens of molars [l^, [4(| . The question 
of whether these effects are significant in formulating phenomenological models 
in cellular physiology, where the typical ionic concentrations are two orders of 
magnitude lower, is largely unexplored. This exploration is beyond the scope of 
the present paper but we point out that our formalism allows the incorporation 



of such effects 41 



The transmembrane flux jf. is now a function of the membrane potential [(/)] 
in addition to the dependencies discussed in the previous section. We require 
that jk satisfy condition (|2.6|) . 



11 



The electrostatic potential <f> satisfies the Poisson equation: 



N 

-V-(eV(f>)=Y / qz k c k (3.2) 

k=l 

where e is the dielectric constant. We shall assume that e is constant in space and 
time. This restriction may be lifted, at the expense of introducing a relation that 
describes the evolution of e. We also assume that there is no fixed background 
charge. It is easy to generalize the calculations below to the case when the 
immobile charges, if present, always stay away from the moving membrane. 
Otherwise, one would need to introduce "collision rules" to determine what 
happens when the membrane hits the immobile charges. We impose Neumann 
boundary conditions for (|3.2[) on the outer boundary r out for simplicity. On the 
membrane T, we impose the following boundary condition: 



dn 



d0 

r dn 



Cm[# (3.3) 



where C m is the capacitance per unit area of membrane. The above is simply a 
statement about the continuity of the electric flux density. Since the membrane 
is moving, the capacitance C m is itself an evolving quantity. We assume the 
following family of constitutive laws for C rn . At x = X(0, t), 

C m (x) = C m (Q(X)) (3.4) 

where Q(X) is the Jacobian or metric determinant of the configuration T t at 
time t with respect to the reference configuration r ro f. This factor describes 
the extent to which the membrane is stretched from the rest configuration. A 
simple example of (|3.4j) would be: 



C m (x) = C° m = const. (3.5) 

As another example, we may set: 

C m (x) =C*Q(X) (3.6) 

where = const is the capacitance per unit area measured in the reference 
configuration. Relation (|3.6[) is the natural scaling if we assume that the mem- 
brane is made of an incompressible material. For suppose the membrane is made 
of a material whose dielectric constant is e m . If the thickness of the membrane 
at the point x = X(0,i) is c£(x), the membrane capacitance there is given by 
e m /<f(x). The incompressibility of the material implies that the local membrane 
volume remains constant in time: d(x)Q(X) = const. Thus, C m (x) must be 
proportional to Q(X). 

Force balance must be modified to take into account electrostatic forces. The 
flow field u satisfies the Stokes equation in fi, or fi e with an electrostatic force 
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N 




V. 

\k=l 

Note that the above equations can also be written as follows: 

V • (S m (u,p) + E e (0)) = 0, V • u = 0, 

E m (u,p) = v{Vu + (V«) T ) - pi, 

S e (0) = e[v0®V0-i|V0| 2 / 



(3.7) 



(3.8) 



Here, E e is the Maxwell stress tensor generated by the electric field. Note that 
we have used (|3.2j) to rewrite the electrostatic force in (|3.7[) in terms of E e . 

We now turn to boundary conditions. Wc continue to let u = on the outer 
boundary r out . On the cell membrane T, we have the following conditions. 
First, by force balance, we have: 

(u,p) + £ e (0))n] = F cl 

as ~t~ F ca p (3.9) 

In addition to F c i as , wc have an additional term F cap which arises because the 
membrane carries capacitive energy. We shall call this the capacitive force, 
which is given as: 



Fcap = T cap K r n - V r T cap , T cap = ^ (c m + ) 



where Kr is the sum of the principal curvatures of the membrane T and Vr = 
V — n(n • V) is the surface gradient on T. The above expression shows that the 
capacitive force can be seen as a surface tension of strength — r cap . The above 
capacitive force is chosen so that Theorem 13 . 1 1 holds . and in this sense, the proof 
of Theorem 13. II provides a variational interpretation of this force. In Appendix 
|Appcndix A.l| we give a physical interpretation of expression (|3.10p . 

An interesting variant of (|3.9[) is the following. Suppose the membrane 
is incompressible in the sense that Q = 1 for all time. We note that this 
condition of two-dimensional incompressibility is not the same as assuming that 
the membrane is made of a (three-dimensional) incompressible material. In the 
case of three-dimensional incompressibility, the membrane may stretch, but this 
would lead to a thinning of the membrane, leading to the constitutive law (|3.6[) 
as we saw earlier. When Q = 1 for all time, we let: 

[(£ m (u,p) + £eW>))n] = Fdas + F p (3.11) 

where F p is given as: 

F p = A/ipn — VpA. (3.12) 

The above is a surface pressure and A is determined so that Q = 1. Note that 
in (|3.11[) we do not need a capacitive force since it can be absorbed into the 
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surface pressure term. 

The continuity condition (|2 . 14[) remains the same. We continue to require 
that the passive flux jk satisfy (|2.6|) . An important difference, however, is that 
jk now depends strongly on the membrane potential [<j>] , given that fXk depends 
on <fi. Ions usually flow through ionic channels, which often have open and closed 
states. The passive flux jk may also depend on such states, which are described 
by gating variables. This flux is the subject of the large experimental and 
theoretical work on ion channels [l2j . For the present work, it is appropriate to 
use the classical phenomenological treatments of flux, although their molecular 
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underpinnings are not clear [42J, l43l |44| . Some of popular choices jh include 



J™ - 9k[»k] = 9k + In l^^j J , (3.13) 

4 GHK = W ( c ^M^') -c k \ T { l = M, (3.14) 
Jk ^ exp(z fe 0')-i ) ,V k B T V ; 

where gk and Pk are positive and depend on the gating variables in certain 
modeling contexts. It is easily seen that both j™ and jJ^ HK satisfy (|2.6[) . We 
shall use expression (|3.14[) in our numerical computations in Section [5l 

We remark that the model we just proposed is nothing other than the 
Poisson-Nernst-Planck-Stokes system if we let uj = ujq given in (|2.ip (iij . The 
novelty here is in the interface conditions at the membrane, (|2.5p . (|3.3p , p.9p 
and (|2.14p . The Poisson Nernst-Planck system has received much attention in 



the field of semiconductors [47], |48|, |49[ , ionic channels [5fJ ,_ion exchange mem 



branes and desalination [46[ as well as physical chemistry [51| 

3.2. Free Energy Identity 

We now show that the system described in the previous section possesses a 
natural free energy. 



Theorem 3.1. Suppose Ck,u,p and (j) are smooth functions that satisfy 
(|3.7p . and (|3.2p in fli andQ e and satisfy boundary conditions (|2.5p . Q3.9p . (|2.14p 

and p.3p on the membrane T. Suppose further that Ck and <p satisfy no-flux 
boundary conditions and u = on the outer boundary T out . Then, Cfc,u,p and 
4> satisfy the following free energy identity. 

"77 (Gs + Eelas + E e i ec ) = —I p — J p — J a 



dt 



E elec = j \e\V<t>\ 2 dx+ f \c m [$\ 2 dm T 



Here, G,s, E e i aSl I p , J Pl I a are the same as in (|2.20p . The same identity holds if 
we require Q = 1 and adopt ()3.1ip instead of p.9p . 
If ak = 0, the free energy is monotone decreasing. 
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In addition to the terms present in Q3.15p . we now have an electrostatic term 
in the energy. 

Before proving Thcorcm l3.1[ we collect some calculus results. We first intro- 
duce some notation. Take a point x = xo G T at t = to- Let X n (t; x$, to) be the 
space-time curve that goes through x = xo at time t = to and is orthogonal to T 
at each time instant. Equivalcntly, X n (i;xo,io) is the solution to the following 
ordinary differential equation: 

^-X n (t; x , t ) = v r (X n , t)n(X n , t), X n (t ; , x , t ) = x . (3.16) 
at 

Here, n(x, t) is the unit normal at the point x at time t pointing from fli into 
fl e , and «r(x, t)n(x, t) is the normal velocity of T at that point. Consider a 
quantity w(x, t) defined on the evolving surface V. Define: 



(U»(xo,to)= ^ W (X n (t;x ,t ),t) 



(3.17) 

x=x ,t=to 



The above expression is an analogue of the convective derivative on the surface 
T. We shall make use of the following well-known identity: 



— / wdmr 
dt 1 



J (D?w + n r wvr) dm r (3.18) 



where kt is the sum of the principal curvatures of T. We now state two calculus 
identities that we shall find useful in the proof of Theorem 13.11 



Lemma 3.2. Let u>(x, t) be a smooth function on T t . We have: 

fr ( wQl ljt) dmr = L (KrU,n ~ (Vrw)) ' ^dT dmr - (3 ' 19) 

where Q is the Jacobian determinant of T t with respect to the reference config- 
uration T re f. 

Proof. Note that 

— = Dfw + {Vrw)-— (3.20) 

where the partial derivatives in t is along material trajectories (constant 0). 
The validity of the above identity should be clear by considering the geometric 
relation between the orthogonal trajectory X n and the material trajectory X. 
We also have the following relation for the time derivative of the integral of w 
over r. 



— [ wdmr = — 
dt r dt 



'dw n -idQ , 



[ wQdm Tlc{ = [ ( ijtQ + ) dmr iet 
k ei Jr,Adt dt) (32i) 
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Comparing this with (|3.18p (with Vr = • n) and using the identity (|3.20[) , we 
obtain the desired result. □ 

Lemma 3.3. Suppose w(x, i),x £ (O, U T) is a smooth function defined in 
Cli whose derivatives are continuous up to the boundary T. Then, we have the 
following identity: 



WD * (In) + [ KvW ^n ~ |Vrw|2 



(3.22) 



where J r denotes integration over the fli face ofT. A similar identity holds for 
functions defined in Q e U T. 

As we shall see, we only need w to be defined in the vicinity of T for the 
above to be true. 

Proof. We only treat the I\ case. The proof for T e is exactly the same. We 
decompose the integrand in the left hand side of (|3.22[) into tangential and 
normal contributions. It is well-known that the Laplacian can be written as: 

. d 2 w dw . .„ „„. 

Aw = — -^ + K r — + A r w (3.23) 
an^ an 

where Ar is the Laplace-Beltrami operator of the surface T. 

We now rewrite §- t (|^) in (l3~22l) in an analogous fashion. For this, we first 
introduce the signed distance function ipix, t) in a neighborhood of T: 

'dist(x,r t ) if x 6 O e , 

{0 ifxer t , (3.24) 

-dist(x, r t ) ifxeO,, 

where dist(x, T t ) is the distance between x and T t - Clearly, Vip evaluated at 
any point on T gives the outward unit normal vector n. Introduce the following 
vector field v defined in a neighborhood of T where ip is smooth: 

v = v r n on T, (Vv)V?/> = 0. (3.25) 

The second condition above just says that v is constant along lines perpendicular 
to the level sets. It is well known that the signed distance function satisfies the 
following transport equation in a neighborhood of T: 

D v ip = ^r- + v • Vtp = 0. (3.26) 
at 

Note that the above convective derivative evaluated on V is equal to Df defined 
in (|3T71) . 
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For any point on T: 

A(^)=V*.V„« (3,7) 

where the subscript t indicates the partial derivative with respect to t. We now 
rewrite this expression as follows: 

W • Viy t = D v (Vip • Vw) - VV>t • Vw - v • V(VV> • Vw) 

= D v (Vip ■ Vra) + V(v • W) • Vw - v • V(V-0 • Vw) 

where we used (|3.26p in the last equality. Now, consider the second term in the 
last line: 

V(v • Vip) • Vw =V r (v • Vip) • V r w + (Vip • V(v • Vip))(Vtp • Vw) 

=V r (v • Vip) • V r w + (V-0 • ({Vv)Vip))(Vip • Vw) (3.29) 
+ (v((V 2 ip)Vip))(Vip-Vw) 

where Vr is the surface gradient on T. Note that V 2 ip is «oi the Laplacian but 
the matrix of second derivatives of ip. The second to last term in (|3.29p is by 
(|3.25p . The last term is also 0, since: 

(V 2 ^)W=iv(|V^| 2 ) = (3.30) 

where we used \Vip\ 2 = 1. Thus (|3.29j) reduces to 

V(v- Vip) -Vw = V r (v- Vip) ■ Vrw. (3.31) 
Let us look at the final term in (I3.28|) . 

v-V(Vf Vw) = (v-Vip)(Vtp-(Vip-Vw)) = (v-Vip)(Vtp-((V 2 w)Vip)) (3.32) 

where we used (|3.30p in the last equality. Combining (|3.3ip . (|3.32p and (|3.28|) . 

we have: 

Vtp-Vw t = D v (Vip-Vw) + V r (v-Vip)-V r w-(vViJj)(Vip-((V 2 w)Vip)) (3.33) 
or equivalently: 

d ( dw\ / dw\ _ „ f d 2 w\ ,„„,•> 

toUJ=^(fa) + Vrwr • Vrw - wr M (3 - 34) 

where we used (|3.25j) . n = Vip on T and the equality of Df and on V. 
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Now, consider the integral: 

J Tt [ W -L (^) + ( wAw ^) dm ? 
= J (wDf ^ + u^Vrwr ■ Vru> + w (atw + K r^— ^ wr^j cfoir (3.35) 

where we used (|3.23p and p. 341) in the first equality and integrated by parts 
along r in the second equality. Note that there are no boundary terms since T 
is a closed compact surface. This proves (|3.22[) . □ 



Wc are now ready to prove Theorem 13.1 



Proof of Theorem \3.1\ First, multiply <\2A\ with fj,k in (|3.1[) and integrate over 
fli and sum in k: 



k— 1 /c— 1 ^ 

The summand in the right hand side becomes: 

MfcV • (ck-^^V^k) dx= [ ( AifeCfe^^V^fe • n ] dm r 

(3.37) 



\ k B T J J T \ k B T 



where n is the outward normal on I\ Consider the left hand side of ()3.36j) . 

Edc k ^ ( dc k i dc k\ 

(3.38) 

fc=i \fc=i / 



We used (|3.1I) in the first equality and (|3.2[) in the last equality. Integrate final 
expression in (|3.38[) over Qj. 



c- 



i r . \ dn \ dt J J J n . dt 



dmj 



g , (3-39) 



(W||V0| 2 ) dx. 
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For the second term in the left hand side of (|3.36[) . we have, similarly to (|3.38|) : 

JV / N \ 

J2 MfcU • Vc ft = V • (ucj) + V I uJ2l z kCk) • (3.40) 

k=l \ k=l / 

Integrate the above expression over f^: 




(3.41) 



Collecting the above calculations, we have rewritten identity (|3.36[) as: 




(3.42) 

where we used qz^cj) = fXf, — (Eq. (|3. 1|> ) in the boundary integral after the 
equality Performing a similar calculation on f2 e , and adding this to the above, 
wc find: 




(3.43) 



where we used (|2.3j) and (|2.5[) to rewrite the boundary integral after the equality. 
Note that the second boundary integral before the equality comes from the fact 
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that the boundary T is moving. Rearranging terms and using (|2.14p . we have: 
f (W^|V0| 2 )dx 



d_ 
~dt 

L 



d 



|V0| 2 + 0V • (eV0) 



dn \ dt 

N \ 

[ir w ]j w + y^[//fe](jfc + a k )\ dm T 



dX 
~dt 



n dm 



r 



(3.44) 



+ / f — c fc - — ^ |V// fc | 2 + (y^qzhCk) u-V0] <ix. 

We used \ik — a k = gzfc</i and used (|3.2[) to rewrite the first boundary integral 
after the equality. Note that: 



8 



dn \ dt 
\ dn 



<9X 

e "H~ ) (eV0)— -n 



Kr^e- e |V r <; 

an 



"si" 



(3.45) 



dmr 



where we used Lemma 13.31 with w = cf> and vr = ^ ■ n in (|3.22j) . Using this 
and the definition of Vr, we may rewrite the first boundary integral in (|3.44|) 
as: 



d ( d<t> 



dn \ dt 

[<j>]D?{C m [(j>])dl 

( -«rC m [0] 5 



||V0| 2 +</>V-(eV0) 



dn 



-|Vr 



— ■ ndm T 
dt 



(3.46) 



where we used (|3.3|) . 

We now turn to equation (|3.7|) . Multiply this by u and integrate over fi^: 

/ u • (cAu — Vp)dx — / ?2fcCfe J u • V</>dx 

= (£(u,p)n) -udmv - / ^|Vu| 2 dx (3.47) 

- / gzfcCfc u • \7(f>dx = 
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Performing a similar calculation on f2 e and by summation, we have: 



[(E(u,p)n) ■ u] dmi 



v IVul dx. 



(3.48) 

Let us first assume (|3.9j) holds. First write S e (0)n in the following form: 

\ (3-49) 



( ^V0-i|V^| 2 n 



00 



We may now write (|3.9p as: 
[£ m (u,p)n] = F olas + F cap 



9n 



90, 



- |V r 0| I n+^Vr( 



9o 



9n 



iv r </>r 



n + C m [0]V r M (3-50) 



where we used (|3.3[) in the last term. 

Using (f3~46]) . (|3T48| and f[3~50]> in (f^44|) we have: 

sL n > + i' v *'> 

= ^ f-[0]-D t n (C7 m [</>]) + (-^ r CU0] 2 n + C m [0]V r [0] + F 
+ y Fdas ■ ^-tor - J \ [i> w ]jw + y^[Mfc](jft + afc)^j d?™r 



ax 

apj 



c?mr 



AT 



(3.51) 



Comparing (|3.5ip , Q3.15P and using (|2.13p . we see that the proof of (|3. 15[) rests 
on the evaluation of the first boundary integral after the equality in (|3.5ip . 



[4>]D?{C m [<f>])dmr = J fof (^C m [0f^j + \{D?C m )[4>A dm T 



1 

dt J r 2 



C m [ci)} 2 dmi 



(3.52) 



where we used (|3.18p with w = [(/>] , = • n in the second equality. We also 
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have: 

(C m [0]Vr[0])dmr = jf (V r Q^mfe] 3 ) - ^(VrC m )[</>] 2 ^ rfm r (3.53) 
Using (|3.52p and (|3.53[) . we have: 

j ^-[^DfiCn^]) + (-KrC m M 2 n + C m [0]VrM) ■ ^ dmr 

= ~ Jt St \ Cm[, ^ dmT ~ S r \ ( D t Cm + (VrCm) ' ^r) [ ^ 2dmr 
+ / (- (^mWj Krn + Vr QcUfl 2 )) ' ^ rf ™r 
Consider the second boundary integral after the equality. 



(3.54) 



S \ ( D ? C ™ + (VrC m ) • ^) [<f\ 2 dm T 
l 9C m ,^ 2 \^ f fldC m 8Q T , 



2 at ^) d ^ = J T {2^i^ ] ) dmT (3 - 55) 

where we used (|3.20[) with u; = C TO in the first equality, and (|3.19[) with w = 
\Q^§^[4>] 2 in the last equality. From ([3351) . ([334]) . (j3~5Tj) and expression (j3~TU)) 
of F cap , we obtain the desired result. 

HQ = 1 and (|3.11l) holds, we may argue as follows. Equation (|3.51|) remains 
valid with F cap replaced by F p . Verification of (|3.15j) rests on the evaluation of 
the first boundary integral in p. 51 j) . Proceeding as in the above, we have: 

J (-[0]£> t "(C m [<fl) + (- Kr C m M 2 n + C m [0]V r [0] +F p ) ■ ^\ dm T 
ni -C m [0] 2 dm r - [ (1^^-uAdmT 



dt J T 2 L ^ J J v \2 dQ dt 

+ St (( A ~ \ Cm ^) KrU ~ Vr ( A ~ \ C ^fj) ■ ^ dm r 

(3.56) 

where we used the (|3.12[) . Since ^ = 0, the second boundary integral after the 
equality is 0. Using p,19[) with w = X — ^C m [4>] 2 and ^ = 0, we see that the 
last boundary integral is also 0. 

In the absence of active currents a^, the free energy is decreasing since jk 
and j w satisfying conditions (|2.6[) and (|2.16|) respectively. □ 
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4. Limiting Systems 



Wc now discuss some limiting cases of the model we introduced in the previ- 
ous section. For this purpose, we shall first make the equations dimcnsionlcss. In 
what follows, the primed symbols denote dimcnsionlcss variables. We introduce 
the following non-dimensionalization of space and time. 

x = Lx', X = LX', t = T D t', T D = D k = D D' k , (4.1) 

where L is the characteristic length scale (for example the size of the domain 
fli) and Do is the characteristic diffusion coefficients of ions. We thus measure 
time with respect to the diffusive time scale of ions. For concentrations and the 
electrostatic potential, we let: 

c k = co4 4> = (4.2) 

For pressure and the membrane elastic force, we let: 

p = cok B Tp', F clas = cok B TF' elas . (4.3) 

For the characteristic fluid and membrane velocity, we turn to relation (|2.14[) . 
Let £ be the characteristic hydraulic permeability of the membrane, which we 
may take as follows: 

dju 



(4.4) 



Then, (^coksT is the characteristic velocity generated by an osmotic gradient 
across the membrane. We thus let: 



u = (c k B Tu'. (4.5) 

With the above dimensionless variables, we may rewrite our system as follows. 
For simplicity, we shall adopt expression (|2.1[) as our definition of the entropic 
part of the free energy lu, so that : 

H k = z k cj)' + \nc k . (4.6) 
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In Qi and Q e , we have: 



84 

at 



f + PeV ■ (u'4) = -V • f k , i' k = -D' k (Vc' k + z k J k V4>'), (4.7a) 



JV 



-V'-(/3 2 V>') = E^ c fe' ( 4Jb ) 



k=l 



7 av = vv+ (X>4J v0', v-u' 



(4.7c) 




where V, V- and A' are the gradient, divergence and Laplace operators in the 
x' coordinate and the dimensionless constants are: 

Pc = D/L > P=-' r d=\l -3—' 7 = — ■ 

In the above, Pe is the Pcclet number which, in this case, measures the ratio 
between the fluid velocity induced by osmotic gradients and the characteristic 
diffusive velocity. The constant j3 measures the ratio between r^, the Debyc 
length and L. The constant 7 is the ratio between the viscosity of water and 
the hydraulic resistance of the membrane. The boundary conditions at the 
membrane interface T become: 



a. 



oij't + 4). («c) 



r, 



(^VV-n)| r . e =eC m [4>'], (4.7f) 

1 dX' , 

7T=3>, (4.7g) 



Pe dt' 



[(S:„(u',p') + /3 2 S' e (0')) n] = + /30F^ ap . (4.7h) 

In equation ()4.7e|) , a is a dimensionless constant given by the ratio of the char- 
acteristic membrane permeability p m and diffusion in the bulk: 



p m _ k B T dj k 



(4.8) 

[**k]=0 



The currents jk and a k are scaled so that jf. = p m coj' k and a k = p m coci' k . In 
621): 

C n = C m C' m ,6= C ° mkBT/q , (4.9) 

qcQr d 

where is the characteristic magnitude of the membrane capacitance per unit 
area (see (|3.5j) or (|3.6j) ). The dimensionless constant 9 is the ratio between the 
membrane charge and the total amount of charge in a layer of thickness on the 
order of the Debye length. In Q4.7g| ), j w = (cokBTj' w . The variables in (|4.7h|) . 
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are defined by: 



£' m (u>') = 7 (V'u' + (V'u') 1 )-p'I, (4.10) 
SU^) = V'0'®V^-i|VV| 2 /, (4.11) 



1 / dC 

2 I m + y dQ 



Fcap = 4 /r n - V^ ap) r^ ap = ^ ( C' m + ) (4.12) 



where k' t (= ktL) is the sum of the two principal curvatures of T measured in 
the x' spatial variable and V r is the surface gradient operator with respect to 
x'. Equations (|4.7a[) - (|4.7c[) and the boundary conditions (|4.7cl) - (|4.71a[) constitute 
our dimcnsionless system. In the rest of this section we shall drop the / in the 
variables with the understanding that all variables, unless otherwise stated, are 
dimcnsionless. 

The dimensionless system above possesses five dimensionless constants a, /3, 7, 
and Pe. We consider two limiting cases of the above system. First of all, con- 
sider the case when Pe <C 1. Assuming all primed quantities are 0(1) with 
respect to Pe, we see, from Q4.7gP that 

— = C(Pe). (4.13) 

Therefore, the membrane does not move to leading order. If we collect all 
leading order terms, we see that the equations (|4.7a[) and (|4.7b[) decouple from 
(|4.7c[) . We thus obtain the following Poisson-Nernst-Planck system with inter- 
face boundary conditions: 

^ = -V • f fc in Q ite (4.14a) 

N 

-V • (/3 2 V(/») = z kC k in n %e (4.14b) 
fe=i 

(f fe -n)| r ._ e =a(j k + a k ), (4.14c) 

-(PV4>'n)\ rtia =0C m [<f>], (4.14d) 

where the membrane T is fixed in time. This model was introduced in [52| (see 
also [53, 54 1 for related models). 



For single cell systems, the Peclet number is about Pe « 10 _1 to 10~ 4 . The 
above may thus be a good approximation to the full system in the Tjj time 
scale. In the context of multicellular systems, however, L may be large and 
Pe can reach unity, as can be seen from expression (|4~7d| of Pe. It should be 
pointed out that there are situations in which the representative fluid velocity is 
not dictated by the osmotic pressure, in which case one should adopt a different 
definition for the Peclet number. For example, if we arc interested in blood cells 
in a flow environment, the ambient hemodynamic flow velocity should be taken 
as the representative velocity. 
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We note that (|4.14p also satisfies a free energy equality. 

Proposition 4.1. Suppose Ck and are smooth functions that satisfy (|4.14l) . 
Then, the following equality holds: 



dt 



E 



elec) 



-F c - J a . 



(4.15) 



In the above, Gs, E e i ec , F c and J a are dimensionless versions of the correspond- 
ing quantities in (|2T2T))) . ([2~2l"jl and (pTT5j) . 



Proof. This follows from a simple calculation. 



□ 



In this sense, system (|4.14|) may be seen as the system associated with the 
energy law (|4.15j) where the mechanical energy and dissipation in (|3.15|) is dis- 
carded. 

Wc next consider the limit when /3 <c 1. This limit is motivated by the fact 
that the Debye length is approximately lnm in typical physiological systems, 
far smaller than the typical length scale of interest. By formally letting /3 — > 
in (|4.7j) we obtain the following system of equations: 



dek 
dt 



f Pcu • Vc fc 

JV 

z k ck 

fc=i 

7Au — Vp 



OfcJ = 



1 <9X 
Pe~dt 



-V • f fe in S\ e 

in Q iy e, 

0, V • u = in fl;, e , 

/ dX 
c^Peu- — 

j w ii, [S m (u,p)n] 



+ ffc 



n 



on r. 



(4.16a) 

(4.16b) 

(4.16c) 
(4.16d) 

(4.16e) 



We have discarded all terms in (|4.7[) that involve j3 and have eliminated the 
boundary condition (|4.7f|) . The most important feature of the above system 
is that we have, in place of the Poisson equation (|4.7b[) . the electroneutrality 
condition (|4.16bp . The electrostatic potential <fi thus evolves so that the elec- 
troneutrality constraint (|4.16bj) is satisfied at each time instant. Although </> 
is thus determined only implicitly through the electroneutrality condition, it is 
possible to obtain a PDE satisfied by </> by taking the derivative of ()4.16bp with 
respect to t and using (|4.16a[) : 



= V • (aV0 + b) 



N 



N 



^zlD k c k , b 



fc=i 



z k D k Vc k . 

fc=l 



(4.17) 



We point out that the electroneutrality condition does not imply that A(f> = 
as may be erroneously inferred from (|4.7b[) . In fact, as (3 — > 0, A</> may remain 
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order 1 with respect to (3 while the right hand side of (|4.7b[) will go to like 
1 . This is a common fallacy in applications of the electroneutral limit. The 
boundary condition for this elliptic equation can be obtained by taking the sum 
in k in boundary condition ()4.16d|) : 



N 



aV(j) + h) n 



i',. 



a-k) 



k=l 



Suppose 



N 



d[4>] 



X z kik > o. 



(4.18) 



(4.19) 



fc=i 



The above inequality states that current flowing out of the cell should increase 
if [(/)] increases, and is thus satisfied by biophysically reasonable expressions for 
jk- This inequality is clearly satisfied if j k are of the form (|3.13[) or (|3.14[) (sec 
also (|2.7p ). Condition (|4.19p is necessary for the boundary value problem (|4.17p 
and (|4.18p to be uniquely solvable (up to an arbitrary constant), assuming a k 
is a given function of x (and t). 

In connection to (|4.17p and (|4. 1 8[) . wc perform the following calculation to 
illuminate the nature of system (|4.16p as it relates to (|4.7p . Suppose <fi satisfies 
(|4.7p . By taking the time derivative of (|4.7b[) with respect to t, we obtain: 



V • = /3 2 V^- +aV0 + b 



dt 



N 



N 



(4.20) 



a = X] Z k D kCk, b = (^ Pez fcCfeU + z k D k Vc k ) . 



k=l 



k=l 



We used ()4.7a[) in deriving the above. At the boundary, we may use (|4.7fP and 
(|4~7e)) to find that: 



/3 2 V 



dt 



aV0 + b • n 



dt 



k=l 



dX 



{C m [<f)}) + X ( ZkCk ~Q~i~ ' 11 + aZ kUk + 0-k) 



(4.21) 



If we formally let /3 in (|420]) and (j42H), wc obtain (jlTf) and |4~T8j) 
respectively. For the above limit to be justified, we must require that ^ and 

remain order 1 with respect to /? as /3 — > 0. It is thus only when the 
evolution of 4> and [</>] is sufficiently slow that we can reliably use system (|4.16p 
as as an approximation to (|4.7p . 

We see from (|4.20p and (|4.21[) that there are two other time scales in the 
system besides the diffusive time scale Td- The first is the Debye time scale, 
(3 2 Td- This is the relaxation time scale of deviations from electroneutrality. 
This Debye time scale is too small to be of physiological interest, and we may 
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safely ignore the /3 2 -gf terms except for the very short initial layer that may 
exist depending on initial conditions. The other time scale, /39Td, which we 
shall call the cable time scale, is the time scale over which the membrane po- 
tential [<j>] can change. In excitable tissue, the ionic currents jk can change 
on a time scale comparable to f38Trj- The interaction of rapid changes in jk 
and the capacitive current /3#J^(C m [0]) term lead to cable effects, including the 
propagation of action potentials, which is essential in describing a wide range 
of electrophysiological behavior. Such phenomena are usually described by the 
cable model, in which the intracellular and extracellular media are treated as 
ohmic resistive media [HI, HH . 

Setting /3#Jj(C m [</>]) term to could thus be problematic for certain applica- 
tions. It is thus of interest to develop a model in which the term (3 2 ^ is ignored 
but P9^(C m [(j)\) is not, while retaining electrodiffusive and osmotic effects con- 
tained in the full model. Such a model (without osmotic effects) is proposed 



in 55[ and was applied in |56[ to a problem in cardiology. A key ingredient in 
the derivation of such a model is an analysis of a boundary layer that forms at 
the membrane interface T. This boundary layer, in physical terms, correspond 
to charge accumulation on both sides of the membrane. We refer the reader to 



571 . |58| for this model and its relationship to conventional cable models. 

An important feature of system (|4.16|) is that it satisfies the following energy 
equality. 

Proposition 4.2. Suppose Ck,4>i u and p are smooth junctions that satisy sys- 
tem (|4.16j) . Then, the following equality holds: 

j t {G s + E elas ) = -I p -J p -J a . (4.22) 

In the above, Gs, E e i as , I p , J p and J a are dimensionless versions of correspond- 
ing quantities in (|2.20|) and (|3.15|l . 

Proof. The proof follows from a calculation similar to the proof of Theorem 
I2~T1 □ 

Thus, system (|4.16|) may be seen as the system associated with the energy 
principle (|4.22[) in which the electrostatic energy in (|3.15[) is discarded. 



5. Numerical Simulation of Animal Cell Volume Control 

In this section, we take the problem of cell volume control to illustrate some 
aspects of the model we introduced above. Cells contain a large number of 
organic molecules that do not leak out through membrane. This results in 
excess intracellular osmotic pressure, which may cause the cell to burst. Cells 
have developed countermeasures to prevent this from happening. 

We shall use the electroneutral system f|4. 16[) to study cell volume control. 
We continue to work with the dimensionless equations. To simplify matters, 
we suppose that the cell membrane T and the outer boundary r ou t = dfl are 
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concentric spheres for all time and that the velocity field u only has a radial 
component. Assuming the boundary condition u = on r out we immediately 
see that u = throughout Qj U fi e . We can thus drop equation ()4.16c|) and set 
u = wherever u appears in system (|4.16[) . Assuming further that Ck and <f> 
are functions only of the (dimensionless) radial coordinate r, we have: 



N 



^z fc c fe =0, (5.1b) 



k=l 



for < r < R and R < r < R ou t where R{t) is the radius of the membrane 
sphere V and i? ut = const, is a the radius of the outer boundary sphere r out . 
The boundary conditions are: 



/•' \ an . ,', ( 5 - lc ) 



at r = 0, 
c k ^f + a(j k + a k ) at r = R±, 



1 dR 

~Pe~dt = Ju " ^ = F ° las at r = R ' (5. Id) 

where r = i?± denote limiting values as r approaches R from above or below. 
Boundary conditions at R — R ou t, will be specified later. The elastic force -F e ias 
can now be viewed as a scalar quantity since the force is only in the radial 
direction. 

We now develop a numerical algorithm to simulate system (|5.1[) and apply 
this to animal cell volume control as a demonstrative example. 

We first discuss the numerical algorithm used to simulate system (|5.ip . Con- 
sider (|5.ip in the region a < r < b. First, suppose b < R or a > R. Then, we 
have: ^ 

j t j T 2 c k dr = a 2 f k {a) - b 2 .f k (b). (5.2) 

If we let b = R(t) in the above, we must account for the fact that R(t) is 
changing in time. Using (|5.1aj) and (|5.1cj) . we have: 



d 
dt 



R(t) 

r 2 c k dr = a 2 f k (a) - R 2 (t)a(j k + a k ). (5.3) 



A similar expression is true when a = R(t). The above conservation relations 
will be the basis for our discretization. 

Let At be the time step, and let R n be the position of the membrane at 
t = nAt. We divide < r < R n and R n < r < R ou t into N v equal segments. 
Let 

{ kRn if n < / < /v 

' „n ^ " " (5.4) 

R n + Rou }f v R , if N v + 1 < I < 2N V . v ; 
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The 2-th segment is given by rj_i < r < rj. Of the 2iV„ segments, segments 
1 < I < N v are in the interior of the cell, whereas the the rest are in the exterior 
of the cell. In each segment, we have the concentrations c£ l and the electrostatic 
potential </>™. 

Suppose we are to advance from time (n — l)At to nAt. We use a splitting 
scheme. Each time step is divided into two substeps. In the first substep, we 
advance membrane position: 

R n = R n-1 _ Pe] n-l At ( 5 5 ) 

In evaluating j w , we need the osmotic pressure as well as the elastic force -F e ias, 
both of which are evaluated using at time (n — I) At. For concentrations of ions 
at the intracellular and extracellular sides of the membrane, we use c^^y 1 and 

c k~Nv+i respectively. 

In the second substep, we update the concentrations and compute the elec- 
trostatic potential. We use one step of a backward Euler discretization. We 
first describe our discretization for the intracellular region. Define: 

cUr) = h l *^'<* (5.6) 

Suppose first that R n < J?" -1 . For 1 < I < N v - 1, we discretize (g^) to obtain 
an equation for c]J ; : 



l — l 



(5.7) 



+ 4 7 r((rr_ 1 ) 2 /^-i-W l ) 2 /^)At. 
where fj} l is set to for I = and 

r 7l n f c k.l ~ c k,l-l . c k,l + c k,l §kl ~ ^k,l-l\ r , . , . „ 

ft = - D * \ Axi 2 Ax~ )^ovl<l<N v -l, 

(5.8) 

where Axt = R n /N v . Note that the integral in (|5 . T[) can be evaluated exactly 
given expression (|5.6[) . As for segment I = N v , we view the endpoint = R n 
as having evolved from J2™ _1 , and thus discretize (|5.3[) . We have: 

Att rR"' 1 
^ " " ^3 /on\3\n _ / A„„2„n-1 



r N v — l 

+ ^ ((r^_i) 2 /^-i " (^") 2 "(^ + Hi) At. 

(5.9) 

The important point here is that the upper end point of the above integral 
is i?" -1 and not R n . The total membrane fluxes (R n ) 2 a^ and (R n ) 2 j% are 
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evaluated at time nAt, and are thus functions of c£ N ,c^ N +1 and [<f)] n = 

If R n > R n 1 , the discretized equations are the same as (|5.7p and (|5.9[) 
except that in (j5.9p the upper endpoint of the integral in is R n instead of J?" -1 . 
The fact that the endpoint of the integral is time-dependent in (|5.3p is taken 
into account by the extension of c^~ 1 (r) when r > i? n_1 (see (|5.Gp V 

The final equation we impose is that clcctroncutrality be satisfied in each 
segment: 

N 

ZkC%i = for alU. (5.10) 

k=l 

For the extracellular segments N v + 1 < I < 2N V , we essentially use the 
same discretization as in the intracellular segments. The only difference is in 
treating boundary conditions at the / = 2N V segment. We impose either no- flux 
or Dirichlet boundary conditions. For no-flux boundary conditions, we simply 
let f£ N =0 in (|5.7p for I = 2N V . Suppose the Dirichlet boundary conditions 
are given by: 

Ck(Rout,t) = Cfc, e - (5.11) 

In this case, we set: 

3 1 

c fc,e = °k,e, C%. e — -C^2N V ~ c k,2N v - 1 ■ (5-12) 

For either boundary condition, the electrostatic potential is determined only up 
to an additive constant, and we thus set </>™ = i(j>2 Nv /2 — (j^N^-i/^ = 0- 

For the second substep, we thus have equations (|5.7p . (|5.9p and (|5.10p with 
suitable boundary conditions at r% N = R ou t , which we must solve for c£ ; and 
<fif. This system of nonlinear algebraic equations is solved using a Newton 
iteration where the Jacobian matrix is computed analytically. In all simulations 
reported here, we obtained convergence to within a relative tolerance of 10 -12 
within less than 4 iterations. In particular, the electroneutrality condition at 
each time step was satisfied at each point to within 6 x 10~ 14 mmol/£ for all 
simulation results shown below. 

Note that the discretization is conservative. For example, we have: 

2JV„ 

E -T ((r " )3 - W-i) 3 ) c M = const ( 5 - 13 ) 
i=i 6 

so long as we impose the no- flux boundary condition at r = R out . We have 
checked this property numerically, we achieve conservation of ions to 14 to 15 
digits. This property is very important in studying long time behavior. 

We would also like to comment on our use of the backward Euler scheme 
and the Newton iteration in the second substep of each time step. Rather than 
use a backward Euler step, we may split the second substep further into two 
substeps. In the first substep, one compute the updates of given values of 
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at time (n — 1) At and in the second substep, we update Ck using the updated (f>. 
A variant of this scheme is to use the above as one step of a fixed point iteration 
to solve the backward Eulcr problem. An advantage of these schemes is that 
the associated matrix problem is much simpler and smaller than that of a full 
Newton iteration we use in this paper. This was indeed the first algorithm we 
used in our attempt to simulate the system. This algorithm, however, turned 
out to have serious stability and convergence issues and led to large pile-up of 
charges close to the membrane. This difficulty was clearly caused by the moving 
membrane. Indeed, a similar algorithm was successfully used in [59j to simulate 
a similar but higher dimensional system, in which the membrane was stationary. 
We also found that if At or the membrane velocity is very small, the fixed-point 
algorithm does produce computational results in agreement with those obtained 
using a Newton iteration. We do point out that even the backward Eulcr, 
Newton scheme, that we use here was not unconditionally stable, though the 
time step restriction was never serious. A more stable algorithm may be possible 
by developing a scheme in which the membrane position and concentrations (and 
electrostatic potential) are computed simultaneously. 

We now describe the model example we simulate. The cell membrane of 
animal cells is not mechanically strong enough to resist osmotic pressure due to 
the presence of organic solutes in the cell. Cell volume control is achieved by 
actively maintaining a concentration gradient of ions across the cell membrane. 
Many modeling studies have been performed to study cell volume control in 
animal cells. To the best of our knowledge, all such studies use ODE systems 
in which the cellular and extracellular concentrations are assumed to have no 
spatial variation [13, EH, 0, EH, |62| . The novelty here is that we use the PDE 
system (|5.1j) . a field theory, to study cell volume control. 

We consider a generic spherical animal cell whose sodium and potassium 
concentration differences across the membrane is maintained by the presence 
of the Na-K ATPasc. Henceforth, we shall use variables with their original 
dimensions, since we will be dealing with a concrete biophysical setup. We 
consider four species of ion, Na + , K + , Cl~ and the organic anions, which we 
index as k = 1, • ■ • , 4 in this order. The diffusion coefficients of the four species 
is given in the Tabled We make the simplification that the organic anions are a 
homogeneous species with a single diffusion coefficient. The diffusion coefficient 
for the organic anion is somewhat arbitrary, one order of magnitude smaller 
than the small inorganic ions. 

We take the initial radius of the spherical cell to be Rq. We let the outer 
edge of the simulation domain i? ut = 2i?o- We assume that the membrane 
does not generate any mechanical force, so that -Foias = 0. Passive water flow 
across the membrane is proportional to the water chemical potential. Given 
that -Fdas = 0, water flow across the cell membrane is driven osmotic pressure 
difference across the membrane: 



4 




(5.14) 
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where Na is the Avogadro constant (so that Na^b is the ideal gas constant), 
T is the absolute temperature, c k is measured in mmol/l and C is measured in 
velocity per pressure. 

For the passive membrane flux j k , we take expression (|3 . 14[) : 

go p ,,( c k \ R _eny{z k <P') - c k \ R+ \ q [ft 

R(ty \ exp{z k (f>') - 1 / k B T 

where the subscript R— and R+ denote evaluation at the inner and outer faces 



of the membrane. This choice is standard for cell volume studies [63l. |64| . The 
number P k is measured in cm/sec and is the permeability of a unit area of mem- 
brane for ionic species k when the radius of the cell is Rq. Assuming that this 
permeability is determined by the presence of ionic channels and that the num- 
ber of ionic channels remains constant, j k must be made inversely proportional 
to the membrane area. For sodium, potassium and chloride, P k is positive but 
we set the permeability for organic solutes to 0. 



We follow [64| to use the following expression for the Na-K ATPasc flux: 

«! = A p ( ^ V ( C2 '*+ V, a 2 = -\ ai . (5.16) 
\cx\ R _+K Na ) \c 2 \ r+ +KkJ 3 

Recall here that oi , c\ are the active Na + flux and concentration respectively 
and ai , c 2 are the active K + flux and concentration respectively. The exponents 
of 3, 2 and the factor of —2/3 reflects the 3 : 2 stoichiometry of the NaK ATPase 
in pumping Na + out and K + into the cell. The constants K^ a and K k are given 
in Tableffl 

All constants and initial conditions are given in Tables [T] and [5] Initial 
concentrations arc assumed spatially uniform. The constants that are not listed 
in the tables are computed so that the initial state is a stationary state under 
no-flux boundary conditions at R = i? ut- This is similar to what is done 
in [(Ml ]. This procedure determines the initial intracellular CI - concentration, 
Na-K ATPase maximal pump rate A p , K + permeability p 2 , initial intracellular 
organic solute concentration, and the organic solute valence Z4. We point out 
that [0] mlt , the initial value of the membrane potential is only needed to compute 
the initial conditions. Once all the concentrations are known, the concentrations 
serve as the initial conditions and there is no need to know [<fi] at the initial time 
to evolve the system forward. 

In the simulations to follow, we took N v = 100 and the time step At = 
500ms. 

We perform the following numerical experiments. Starting with the initial 
conditions specified above with no-flux boundary conditions, we set the following 
Dirichlet boundary conditions for t > 10s: 

ci, e = 100, c 2 , e = 50, c 3 , e - 150, (5.17) 

where the units are in mmol/£. The boundary concentrations are thus isotonic 
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T(K) 


273.15 + 37 


Xfe (mmol/f) 


0.75 [64] 


£ (cm/s/mPa) 


5.2507 x 10~ 13 [63] 


KjVo (mmol/^) 


3.5 [64] 


i?o (mm) 


0.5 


(cm/s) 




R out (mm) 


1 


[0] init (mV) 


-70 



Table 1: Constants used in the numerical simulation. [</>] lnlt is the initial membrane voltage. 
Symbols labeled with '-' are determined so that the initial condition is a stationary state (see 
main text). The ion related constants are listed in Table [2] 





Zk 


D k (cm 2 /s) 


P k (cm/s) 


init 
L /c,int 


init 

L fc,cxt 


Na+ 


+1 


1.33 x lO" 5 ^ 
1.96 x 10~ 5 [65] 
2.03 x 10~ 5 [65] 


1.0 x 10- y [66] 


10 


145 


K+ 


+1 




140 


5 


cr 


-1 


1.0 x 10~ 7 [66] 




150 


O.A. 




1.0 x 10" 6 











Table 2: Parameters related to ionic concentrations. O.A. stands for organic anions. The 
initial intracellular and extracellular concentrations are given by cj?j* t and c 1 ™* t respectively 
(listed here in mmol/Q. Symbols labeled with '-' are determined so that the initial condition 
is a stationary state (see main text). Other parameters are listed in Tabic Q] 

with the initial concentrations, but the extracellular K + concentration is now 
increased 10-fold. Such a stimulus should lead to immediate depolarization 
together with expansion of the cell. 

The computational results are given in Figure [TJ What is interesting here is 
that there is a transient drop in the cell radius, followed by an expected gradual 
increase. This transient drop is due to the following. After a sudden change in 
the boundary condition, Na + ions diffuse out whereas K + ions should diffuse in 
from R = i? ut- Since K + diffuses faster than Na + , there is a transient increase 
in total ionic concentration near the membrane, leading to excess osmotic pres- 
sure immediately outside the cell compared to the inside. This gives rise to a 
transient drop in the cell radius. However, as the ionic concentration becomes 
spatially uniform within the extracellular and intracellular domains, the cell 
starts to expand. 

The next computational results describe a hypotonic shock. We set the 
boundary conditions to the following for t > 10s: 

ci, e = 100, c 2 , e = 5, c 3 , e = 105, (5.18) 

where the concentrations are in mmol/£. A snapshot of the computational 
results are given in Figure [2J The cell expands due to the hypotonic shock but 
tends to a new stationary state with time. 
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Figure 1: Computational results under a high K + stimulus (see 1 15. 17H ). The first five figures 
are the snapshots of the ionic concentrations and the electrostatic potential at t = 50s. The 
horizontal axis represents the radius r. The last figure plots the cell radius R(t) as a function 
of time. 
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Figure 2: Computational results under hypotonic stimulus (see (15. 18tl ) . The first five figures 
are the snapshots of the ionic concentrations and the electrostatic potential at t = 100s. The 
horizontal axis represents the radius r. The last figure plots the cell radius R(t) as a function 
of time. 
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6. Conclusion 



We introduced a PDE system of clcctrodiffusion and osmotic water flow 
in the presence of dcformablc capacitance-carrying membranes. The salient 
feature of the model is that it satisfies an energy equality, and thus possesses 
a natural thermodynamic structure. We discussed simplifications of the model 
and applied the electroneutral limit to the problem of cell volume control. 

In the proof of Theorem 12.11 we showed that the van t'Hoff expression for 
osmotic pressure arises naturally, simply through an integration by parts argu- 
ment. This observation seems to be new. It is interesting that, in expression 
(|2.20[) . the mechanical pressure p and osmotic pressure ir w only appear in the 
combination ip w = p + tt w . This is consistent with experimental results in- 
dicating that the effect of osmotic pressure on transmembrane water flow is 
indistinguishable from that of mechanical pressure 67 1. 

The models introduced here are sharp interface models in the sense that the 
membrane is treated as a surface without thickness and the physical quanti- 
ties of interest are allowed to have discontinuities across T. This is in contrast 
to diffuse interface models in which the membrane has some small but finite 
thickness and the physical quantities transition rapidly but smoothly across the 
interface. It should be possible to obtain at least parts of the model by taking 
the thin interface limit of an appropriate diffuse interface (or finite thickness) 
model. This may lead to a simpler verification of the energy identity of The- 
orem 13.11 Establishing such a connection may also help in understanding the 
physical nature of the capacitive force (|3.10j) . The calculations in Appendix 
|Appcndix A.l| may be seen as an initial step in establishing this relationship. 

Given the natural thermodynamic structure of the problem, it is almost 
certainly the case that our model has a variational structure. A variational 
principle for dissipative systems dates back to [H^- This procedure has been 
used successfully in deriving dynamic equations for soft matter systems 6^, 7^1 ■ 
In [4l|, [7l[, a model for non-ideal electrolyte solutions is derived by combining, 
in the spirit of Rayleigh (see [72j ) , the principle of least action with the above 
variational principle for dissipative system. 

The energy identities introduced here provide a natural apriori estimate for 
our PDE system. For related systems, the corresponding free e nerg y identity 
has been used successfully to prove stability of steady states [73, 74 1. It would 
be interesting to see if similar analytical results can be obtained for the system 
proposed here. 

We hope that our model has wide-ranging applications in cellular physiology. 
In principle, our model is applicable to most problems of classical physiology 
[J- A s we saw m Section IH our model admits simplifications when certain 
dimensionless parameters are small. In the short time scale when water move- 
ment is not significant, the system is reduced to the Poisson-Nernst-Planck 
model with interface boundary conditions. This and related models have been 
successfully applied in [s2> [lH, [56| • If the physiological processes of interest are 
slow and happen over a long time scale, the electroneutral limit may be taken. 
This was applied to the problem of cell volume control in Section[S]of this paper. 
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Any serious application of our model will require the development of an 
efficient numerical algorithm. The electrodiffusive part of the problem with 
stationary membranes (without fluid flow) has been treated successfully in [59l 
|76| in a two-dimensional setting. In the model presented here, the membrane 
interface is dynamic. We must therefore solve an electrodiffusive problem in 
a domain with a moving interface across which physical quantities experience 
discontinuities. We have presented successful computations in one-dimension 
for the electroneutral limit in Section [SJ but simulations are bound to be more 
challenging in higher dimension. If a regular mesh is to be used, immersed 
boundary or immersed interface schemes could be a major component of the 
algorithm (H |zzl . 

Many physiological phenomena in which both electrodiffusion and osmosis 
play an important role take place over spatial scales of whole tissues or organs 
rather than the cellular spatial scale we focused on in this paper. Such systems 
include ocular fluid circulation, electrolyte regulation in the kidney or brain 
ionic homeostasis. For such systems, it is important to develop an appropriate 
homogenized model. In the context of cable models, this is known as the bido- 
main model, and has fou nd g reat utility in many contexts, especially in cardiac 
electrophysiology [2! HAS HI- We shall report on a such a multidomain 
model in a future publication. 



Appendix A. Appendix 

Appendix A.l. Physical Interpretation of the Capacitive Force 

Let the membrane be made of an incompressible material. In this case, we 
argued that C m should satisfy (|3.6p . We shall show that the incompressibility 
of the material implies r cap = — C m [cj)} 2 , in agreement with p,10[) . To this end, 
we take a membrane of finite thickness d and consider the limit as the thickness 
tends to 0. Let T be the midplane of this membrane of finite thickness. The 
membrane thus coincides with r as d — > 0. 

Take a point x G T and let n be the unit normal and d the thickness of the 
membrane at x. The stress inside the membrane is given by: 

^mcm £unem ^mcinj 

where £™ om is the Maxwell stress, E™°™ the elastic stress p mem J is the isotropic 
pressure term that enforces incompressibility of the material. We have made 
the simplification that the material can only generate isotropic stresses. 

Let us now consider the limiting behavior of this stress when d is very small. 
To leading order in d, the Maxwell stress tensor inside the membrane is given 
by: 

=e m ^(n®n-±A=C m ^(n®n- \l\ (A.2) 

where e m is the dielectric constant of the membrane. We assumed that the 
electric field inside the membrane is given by [4>]n/d, given that the membrane 
is very thin. We used e m /d = C m in the second equality. 
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Now, let us consider stress balance at x+ |n, the point where the membrane 
touches fl e . Here, we have the following stress balance condition: 

E mcm n = E o Cn ( A 3^ 

where S fie is the stress in fi e . Using (jA.lj) and (|A.2|) S mom n, to leading order 
in d, can be written as: 

S mcm n = (c m ^ - p mcm ^j n. (A.4) 

As d — > 0, E^n must remain finite if there is a finite distinguished limit. 
Therefore, S mem n must remain order 1 with respect to d. In (|A.4p . the term 

C m 24 grows like l/d as d — > 0. The elastic stress stays order 1 with respect to 
d. Therefore, p mem n, must satisfy: 

pmom = Cm ^M (A - 5) 

to leading order in d. 

Take any unit vector t tangent to L at x. We have: 

£ mcm t = -ic m [0] 2 t, (A.6) 
d 

to leading order in d. Multiplying the above by the thickness d of the membrane, 
and taking the limit as d —> 0, we conclude that a "surface tension" of magnitude 
— C m {4>] 2 is generated at the membrane. 

The above derivation suggests the following physical interpretation of ex- 
pression (j3~10l) . The term \C m \$f comes directly from the Maxwell stress. 
More simply put, this tension comes from large coulomb forces squeezing the 
thin membrane. This force must be counter balanced to maintain the mechan- 
ical integrity of the membrane, which is given by an isotropic pressure in the 
case of incompressible materials. This contributes the term \Q^§Q^[ ( f\ 2 to the 
capacitive force. 

Appendix A. 2. Relation to an Immersed Boundary Model of Osmosis 



In this appendix, we explore a model of osmosis proposed in [19l . |24| and 
compare this with the model presented in this paper. 

Consider the system of equations introduced in Section [2] We shall take 
<jj = cj given in (|2.1[) . Let there be only one solute species and suppose that it 
is impermeable to the membrane. Let c be the concentration of this solute. We 
have: 

g + V-(uc) = V-OD(Vc)), (A.7a) 
yAu - Vp = 0, (A.7b) 
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and at the boundary, we have: 

[E(u,p)] = F cias , uc-D^ = c ^-n (A.7c) 

and (|2. 14[) . We impose no- flux and no- flow boundary conditions on d£l. We 
know from Theorem 1 2 . 1 1 that the solutions to (|A.7[) satisfy: 



d f 

— / (k B Tc\nc + E eias )dx 

= - / Dc\V\nc\ 2 dx + / (F clas -n + k B T[c})j w dm r . 
Jn Jr 



(Ai 



Dissipation in the free energy will be guaranteed if, for example, j w is a linear 
function of F c i as • n + fcsT[c] of negative slope. 

In [Til [24| , the authors propose the following immersed boundary model of 
osmosis. We argue why this may be considered a regularization of the above 
model. The concentration c satisfies: 

^ + V.(uc) = V.(u(vc + ^V^)). (A.9a) 

Here, tJj v is a "barrier" potential localized near the membrane: 

^(x)=/ <^(x-X(0,i))dmr ref (A.9b) 
Jr rB[ 

<p v (x) =v~ 3 <Pi f^J , (A.9c) 

where ^i(x) is a positive function of compact support. For simplicity, we shall 
take (p(x.) to be radially symmetric and nonzero only if x < 1. The parameter 
77 > is the width of the barrier potential. As 77 — > 0, one would expect the 
barrier to be so high that the solute will be effectively blocked from crossing the 
membrane. 

The fluid velocity field satisfies the Stokes equation in f2\r with an external 
force term that comes from the presence of solutes interacting with the barrier: 

uAu - Vp = cVVv (A.9d) 

At the boundary, we suppose that the velocity field is continuous and that the 
stress tensor H m (u,p) satisfies the following jump condition: 

[E m (u,p)n] = F clas + F so i, (A.9c) 

/ cVy>„(x-X(0,i))dx. (A.9f) 
Jn 

The force F so i should be seen as a reaction force to the cVip term in (IA.9d|) . 



ax 




de 


-J 
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The membrane velocity ^ and u satisfy relation (|2.14p . The following energy 



dt 

identity is satisfied by this system: 
d f 

— / (k B Tc In c + Solas + c4> v ) dx 
di Jq 



i 



J Dc V^lnc+^j dx + J (F e i as + Fsoi) • nj w dm r . 



(A.IO) 



We leave the verification of this identity to the reader. The above energy will 
be monotonically decreasing if j w is a linear function of F c i as + F so i of negative 
slope. 

Let us now consider the limit as r\ — > 0. First, consider the advection- 
diffusion equation satisfied by the solute. The 77 — >■ limit corresponds to the 
barrier potential becoming infinitely thin and high, and thus, this limit should 
lead to ()A.7a[) in f2\r and the impermeable boundary condition (|A.7cj) for the 
solute. 

To examine the limiting behavior of the Stokes equation rewrite (|A.9d[) and 
(|A.9cl) in the following immersed boundary fashion: 

I/An - Vp = cWlpr, + f B ol + feias- (A.ll) 

Here f = f so i or f i as are surface measures supported on T whose action on a 
test function v is given by: 



(f,v}= J F(X)-v(X)dm r (A.12) 



where F = F so i or F c i as . Equation (|A.11I) is thus to be understood in the sense 
of distributions. The interface conditions (|A.9c|) are now expressed in terms of 
singular force fields. This reformulation of interface boundary conditions is a 
key ingredient in the immersed boundary method [jjj]. 

Given that cVip v and f so i are forces of action and reaction, as r\ — > 0, one 
expects cV-0 + f so i to go to in a distributional sense. The right hand side of 
(jA.lip will thus be reduced to f e i as . Equation (|A.11|) will then be equivalent 
imposing the boundary condition (|A.7c[l at L to the Stokes equation satisfied in 

n\r. 

In view of fOji . ([ATO]) and the discussion after these equations, we would 
like to show that F so \ approaches [c]fceT, the van t'Hoff expression of osmotic 
pressure, as r\ — > 0. If this is the case, we can view system (|A.9[) as giving a 
mechanical interpretation of osmotic pressure. In this view, osmotic pressure 
is generated by solute molecules hitting the impermeable membrane. We now 
sketch a boundary layer calculation that lends credence to this claim. Take a 
point on T and let this point be the origin without loss of generality Assume 
r is locally flat. Let L coincide with the x — y plane and let z be the direction 
normal to T the positive z axis to be pointing into f2 e . Assume furthermore, 
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that I S I = const over this flat region. Consider the solute flux for z > 0: 



J, 



9c c dip^ 



dz k B T dz 

Introduce the stretched boundary layer coordinate r\Z = z: 

9c c 9-017 
1 z dZ k B T dZ 

We conclude, to leading order, that 

c(Z) = c(Z = I) exp(-1> v {Z)/k B T) + o(l). 
In original coordinates, we have, for < z < rj, 



c(x', z) = c(x', 77) exp 



o(l) 



(A.13) 



(A.14) 



(A.15) 



(A.16) 



where x' = (x, y) and o(l) denotes terms that tend to as 77 — )• 0. Likewise, for 
—77 < z < 0, we have: 



c(x', z) = c(x', —77) exp 

Let us now compute F so i at x = y = 0. 
_ 9X 



o(l). 



(A.17) 



■ sol 



90 



X|<T) 



cV(/5^(x)dx 



/ cV</?,,(x)(ix + / cV( / 9 r) (x)dx (A, 

7|x|<r7,z>0 ^|x|<r;,z<0 



18) 



(F s + ol + F-J 



9X 



90 



where we used the fact that Lp v is supported in |x| < ?/. Consider the first term: 

__l 9x 



90 



x|<?j,z>0 



x|<?7,Z>0 



cV(/3^(x)dx 
c(x', 77) exp ( 
c(x', 77) exp 



x|<?7,z>0 



k B T 
k B T 



V(^(x)J dx + o(l) 

J ^-ip v {x', z)j dx'dze z + o(l) 

(A.19) 
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where e 2 is the unit coordinate vector in the z direction. Note that: 



<9X 



86 



y'\<v 



¥>r,(y' , z)dy' 



(A.20) 



near the origin. Here, we used the fact that the membrane is flat and that | §jj | 
is constant. Using this, we have: 



sol 



•II 



c(0, rf) exp 



k B T J d 



ipniz) dze z + o(l) 



-c(0, v)k B T + c(0, 0)k B Texp - 



k B T 



o(l). 



Likewise, we have: 



Ko\ = c (°> -v)k B T - c(0,0)fc B Texp - 



^(0) 
k B T 



Thus, as r) — > 0, we have: 

F so i = (c(0,0-) - c(0,0+))fc B Te z 



(A.21) 



+ o(l). (A.22) 



(A.23) 



as desired. 

We emphasize that the above calculation is purely symbolic. We are not 
attempting to prove that the solutions to (|A.9[) approach the solutions to (|A.7J) 
in the 77 — > limit. 
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